 ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ;%  ncl script to draw locations of measurements 
 ;%  written by Kai Zhang, MPI-M 
 ;%  kai.zhang@zmaw.de 
 ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" 
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
  load "$KAILIBROOT/string_number_500.ncl" 

  plotpath = "~/fig/" 
  plotpath = "./" 
  plotname = plotpath+"location" 
  plotform = "eps" 
  ploteps  = plotname+".eps" 
  plotpdf  = plotname+".pdf" 


  colormap = "psgcap"
  icla = 148
  iclb = 106
  iclc = 159 
  icld = 178
  icle = 109

  L_add_square = False
  L_add_text = False   
  L_add_text = True  
  L_add_dot = False   
  L_add_dot = True  
  R_txFontHeightF = 0.012
 
  small_map = True 
  small_map = False   

  lat_filename = "radon_sitelat.dat"
  lon_filename = "radon_sitelon.dat"
  bad_filename = "radon_sitebad.dat" 
  nam_filename = "radon_sitename.dat" 
  prt_filename = "radon_siteprintname.dat" 
  eur_filename = "radon_siteeur.dat"
  chi_filename = "radon_sitechi.dat"
  usa_filename = "radon_siteusa.dat"
  sam_filename = "radon_sitesam.dat"
 
  lat  = asciiread(lat_filename,-1,"float") 
  lon  = asciiread(lon_filename,-1,"float") 
  kbad = asciiread(bad_filename,-1,"integer") 
  ieur = asciiread(eur_filename,-1,"integer") 
  ichi = asciiread(chi_filename,-1,"integer") 
  iusa = asciiread(usa_filename,-1,"integer") 
  isam = asciiread(sam_filename,-1,"integer") 
  kset = kbad
  kset = 1 

  sitename = asciiread(nam_filename,-1,"string") 
  prtname  = asciiread(prt_filename,-1,"string") 
  
printVarSummary(sitename) 
print(prtname) 

  ns = dimsizes(lat) 

  print("") 
  print(" dimension of lat " + dimsizes(lat) ) 
  print(" dimension of lon " + dimsizes(lon) ) 
  print("") 


  wks = gsn_open_wks(plotform,plotname)
  gsn_define_colormap(wks,colormap)
 
  minlat   = -90
  maxlat   =  90.
  minlon   = -180.
  maxlon   =  180.
 
  mpres             = True
  mpres@gsnFrame    = False 
  mpres@gsnMaximize = False 
  mpres@gsnMaximize = True 

  nol = string_number_500(0:dimsizes(lon)-1)  

;-------------------------------------------------------------------------------------------------
; map 
;-------------------------------------------------------------------------------------------------

  mpres@tiMainString = ""
  
  mpres@mpProjection = "CylindricalEquidistant" ;"Mollweide"   

  minlat   = -90
  maxlat   =  90.
  minlon   = -180.
  maxlon   =  180.
   
  mpres@mpMinLatF   = minlat
  mpres@mpMaxLatF   = maxlat
  mpres@mpMinLonF   = minlon
  mpres@mpMaxLonF   = maxlon

  nc1 = NhlNewColor(wks,.8,.8,.8)
  mpres@mpFillColors = (/-1,-1,nc1,-1/) 
 
  peru_rgb  = (/205,133, 64/)/255.   ; RGB triplets for the 
  coral_rgb = (/255,127, 80/)/255.   ; named colors "peru",
  bwood_rgb = (/255,211,155/)/255.   ; "burlywood", and "coral".

  peru   = NhlNewColor(wks,peru_rgb(0),peru_rgb(1),peru_rgb(2))
  coral  = NhlNewColor(wks,coral_rgb(0),coral_rgb(1),coral_rgb(2))
  bwood  = NhlNewColor(wks,bwood_rgb(0),bwood_rgb(1),bwood_rgb(2))
 
  colors   = (/peru,coral,bwood/)         ; fill colors for land
  patterns = (/  2,   7,   14  /)         ; patterns for grid lines

 ;mpres@mpGridAndLimbOn        = True        ; turn on lat/lon lines
 ;mpres@mpGridAndLimbDrawOrder = "PreDraw"   ; draw lines first, so
                                             ; filled land is on top
  mpres@mpPerimOn              = False       ; turn off map perimeter
  mpres@mpOutlineOn            = False       ; turn on map outlines
  mpres@mpGeophysicalLineColor = "Gray21"    ; map outline color

 ;mpres@mpGridLineDashPattern = patterns(1)      ; grid line pattern
  mpres@mpLandFillColor       = "grey" ;colors(2)        ; land fill color
 
 ;mpres@mpGridLineColor       = "Gray21"   
 ;mpres@mpGridLineColor       = "grey" 

;-------------------------------------------------------------------------------------------------
; Tick Mark
;-------------------------------------------------------------------------------------------------

  ;mpres@tmXBMode = "Explicit"         ; explicit labels
  ;mpres@tmXBValues = (/1,2,3,4,5,6,7,8,9,10/)
  ;mpres@tmXBLabels = (/"HU","BE","CH","XI","WU","GU","NA","SH","FU","GA"/)
  
  
  mpres@tmBorderThicknessF = 2.

  ;mpres@tmXBLabelFont = "helvetica"
 
  mpres@tmXBLabelFontHeightF = 0.012
  mpres@tmXTLabelFontHeightF = 0.012
  mpres@tmYLLabelFontHeightF = 0.012
  mpres@tmYRLabelFontHeightF = 0.012
  
  mpres@tmXBLabelFontThicknessF = 0.8
  mpres@tmXTLabelFontThicknessF = 0.8
  mpres@tmYLLabelFontThicknessF = 0.8
  mpres@tmYRLabelFontThicknessF = 0.8

  mpres@tmXBOn = True
  mpres@tmXTOn = True
  mpres@tmYLOn = True
  mpres@tmYROn = True

  ;mpres@tmXBMode = "Manual"
  ;mpres@tmXBTickStartF = 20.
  ;mpres@tmXBTickEndF = 80.
  ;mpres@tmXBTickSpacingF = 5.
  
  ;mpres@gsnMajorLatSpacing = 5.
  ;mpres@gsnMajorLonSpacing = 5.
  ;mpres@gsnMinorLatSpacing = 1.
  ;mpres@gsnMinorLonSpacing = 1.
  
  mpres@tmXBLabelsOn = True
  mpres@tmYLLabelsOn = True

  ;mpres@tmXTLabelsOn = True
  ;mpres@tmYRLabelsOn = True
  
  mpres@tmXBBorderOn = True
  mpres@tmXTBorderOn = True
  mpres@tmYLBorderOn = True
  mpres@tmYRBorderOn = True
  
  mpres@tmXBLabelStride = 1
  mpres@tmXTLabelStride = 1
  mpres@tmYLLabelStride = 1
  mpres@tmYRLabelStride = 1
  
  ;mpres@tmYMajorGridThicknessF = 1.0 
  ;mpres@tmXMajorGridThicknessF = 1.0 

  mpres@tmXBMajorThicknessF = 2. 
  mpres@tmXTMajorThicknessF = 2.
  mpres@tmYLMajorThicknessF = 2.
  mpres@tmYRMajorThicknessF = 2. 
  ;
  mpres@tmXTMajorLengthF = 0.008 ;0.012
  mpres@tmXBMajorLengthF = 0.008 ;0.012
  mpres@tmYLMajorLengthF = 0.008 ;0.012
  mpres@tmYRMajorLengthF = 0.008 ;0.012
  
  mpres@tmXTMinorLengthF = 0.004
  mpres@tmXBMinorLengthF = 0.004
  mpres@tmYLMinorLengthF = 0.004
  mpres@tmYRMinorLengthF = 0.004
  
  mpres@tmXTMajorOutwardLengthF = 0.008 ;0.012
  mpres@tmXBMajorOutwardLengthF = 0.008 ;0.012
  mpres@tmYLMajorOutwardLengthF = 0.008 ;0.012
  mpres@tmYRMajorOutwardLengthF = 0.008 ;0.012
  
  mpres@tmXTMinorOutwardLengthF = 0.004
  mpres@tmXBMinorOutwardLengthF = 0.004
  mpres@tmYLMinorOutwardLengthF = 0.004
  mpres@tmYRMinorOutwardLengthF = 0.004
   

;-------------------------------------------------------------------------------------------------
; shape 
;-------------------------------------------------------------------------------------------------
;   mpres@mpShapeMode  = "FreeAspect"
;  
;   mpres@vpXF = 0.1
;   mpres@vpYF = 0.8
;
;   mpres@vpWidthF = 0.8
;   mpres@vpHeightF = 0.8 * .618

;-------------------------------------------------------------------------------------------------
; plot map 
;-------------------------------------------------------------------------------------------------
  
  map = gsn_csm_map_ce(wks,mpres)
  


;-------------------------------------------------------------------------------------------------
; small map (annotation) 
;-------------------------------------------------------------------------------------------------
 
  if(small_map) then 

  npres                   = True
  npres@gsnFrame          = False
  npres@gsnDraw           = False
 
  npres@pmTickMarkDisplayMode = "Always"
  npres@tmXBLabelFontHeightF  = 0.01
 
  npres@mpDataBaseVersion  ="MediumRes"
  npres@mpOutlineOn       = True
  npres@mpMinLatF         =  40. 
  npres@mpMaxLatF         =  60.
  npres@mpMinLonF         =   0.
  npres@mpMaxLonF         =  20.
  
  npres@vpHeightF         = 0.2    ; Make this second map much smaller.
  npres@vpWidthF          = 0.2

  npres@tmXBOn = False
  npres@tmXTOn = False
  npres@tmYLOn = False
  npres@tmYROn = False

  map2 = gsn_csm_map_ce(wks,npres)

;
; Set up some resources to indicate where the smaller map
; will be (upper left edge of plot).
;
  amres                  = True
  amres@amParallelPosF   = -0.48    ; -0.5 is the left edge of the plot.
  amres@amOrthogonalPosF =  0.03    ; -0.5 is the top edge of the plot.
  amres@amJust           = "TopLeft"

  map_anno = gsn_add_annotation(map, map2, amres)  ; Attach map to map.

  end if 

;-------------------------------------------------------------------------------------------------
; Tick Mark
;-------------------------------------------------------------------------------------------------


  if (L_add_square) then 

  print("  L_add_square = " + L_add_square ) 

  minx = 15.
  maxx = 50.
  miny = 100.
  maxy = 135.
  
  xpts = (/1.,1./)
  ypts = (/1.,1./)
  
;-------------------------------------------------------------
; Add Longitudes
;-------------------------------------------------------------
  resl = True 

  resl@gsLineColor = "black" 
  resl@gsLineThicknessF = 1.0
  resl@gsLineLabelString= ""
  resl@gsLineDashPattern = 5
   
  xlonbeg = miny
  xlonend = maxy
  ylatbeg = minx
  ylatend = maxx

  ypts(0) = ylatbeg
  ypts(1) = ylatbeg 
  xpts(0) = xlonbeg
  xpts(1) = xlonend 

  dlma = gsn_add_polyline(wks,map,xpts(:),ypts(:),resl)      

  ypts(0) = ylatend
  ypts(1) = ylatend 
  xpts(0) = xlonbeg
  xpts(1) = xlonend 

  dlmb = gsn_add_polyline(wks,map,xpts(:),ypts(:),resl)   
  
;-------------------------------------------------------------
; Add Latitudes
;-------------------------------------------------------------
  ypts(0) = ylatbeg
  ypts(1) = ylatend 
  xpts(0) = xlonbeg
  xpts(1) = xlonbeg 
  
  dlmc = gsn_add_polyline(wks,map,xpts(:),ypts(:),resl)      
   
  ypts(0) = ylatbeg
  ypts(1) = ylatend 
  xpts(0) = xlonend
  xpts(1) = xlonend 
  
  dlmd = gsn_add_polyline(wks,map,xpts(:),ypts(:),resl)      
  

;-------------------------------------------------------------
; L_add_square end 
;-------------------------------------------------------------
  end if 


 ;drawNDCGrid(wks)


  ;;print(prtname) 


;-------------------------------------------------------------
; text location 
;-------------------------------------------------------------
  if (L_add_dot) then  
    lonx = where(lon.lt.0,   lon+10.0,lon-10.0) ;;   0-180 left 
    lonx = where(lon.gt.180.,lon+10.0,lon-10.0) ;; 180-360 right 
  else  
    lonx = lon  
  end if

  latx = lat 


;-------------------------------------------------------------
; retune text location 
;-------------------------------------------------------------
  load "plot_radonsite_name.ncl" 


  print(prtname) 
  
  lonx = where(kset.eq.0,lonx,lon) ;; if it is not set yet, use lon 
  latx = where(kset.eq.0,latx,lat) ;; if it is not set yet, use lat+2.0 

;-------------------------------------------------------------
; Add text
;-------------------------------------------------------------


  ;;print(prtname) 
  
  if (L_add_text) then 

  print("  L_add_text = " + L_add_text ) 

  txres = True 
  txres@txFont = "helvetica" ; "times-roman", "helvetica-bold"
  txres@txFontHeightF = R_txFontHeightF

  ik = 0 
  il = 0 

  do is = 0, ns-1  

     if(sitename(is).eq."ispra") then 
        kbad(is) = 0 
     end if 

     if(sitename(is).eq."pathankot") then
        kbad(is) = 0
     end if

     if(sitename(is).eq."new_york") then
        kbad(is) = 0
     end if

     if(sitename(is).eq."toulouse") then
        kbad(is) = 0
     end if

     if(sitename(is).eq."puy") then
        kbad(is) = 0
     end if

     if(sitename(is).eq."jungfraujoch") then
        kbad(is) = 0
     end if

     if(sitename(is).eq."moscow") then
        kbad(is) = 0
     end if

     if(kbad(is).eq.1) then 

        if(kset(is).eq.1) then 
           txres@txFontHeightF = 0.007
           txres@txFontColor = 26
           print(" id, site : " + il + "  " + sitename(is) )  
           ;text = gsn_add_text(wks,map,prtname(is),lonx(is),latx(is),txres)  
           il = il + 1 
        else

       ;print(" id, site : " + ik + "  " + sitename(is) ) 


        ;; smaller the font for some sites 

        txres@txFontHeightF = R_txFontHeightF

        if(False) then 

        if(sitename(is).eq."saclay") then 
        txres@txFontHeightF = R_txFontHeightF * 0.8
        end if 
        if(sitename(is).eq."rovaniemi") then 
        txres@txFontHeightF = R_txFontHeightF * 0.8
        end if 
        if(sitename(is).eq."joensuu") then 
        txres@txFontHeightF = R_txFontHeightF * 0.8
        end if 
        if(sitename(is).eq."nurmijarvi") then 
        txres@txFontHeightF = R_txFontHeightF * 0.8
        end if 
        if(sitename(is).eq."beijing") then 
        txres@txFontHeightF = R_txFontHeightF * 0.8
        end if 

        end if 
        
        tstr = prtname(is) 

        txres@txFontColor = "black" 
       ;text = gsn_add_text(wks,map,nol(ik),lonx(is),latx(is),txres) 
        text = gsn_add_text(wks,map,tstr,lonx(is),latx(is),txres) 
        ik = ik + 1 
        end if 

     end if 
  end do

  draw(map)  
  
  end if  

  
;-------------------------------------------------------------
; Add dot
;-------------------------------------------------------------

  print(iusa) 

  if (L_add_dot) then 

  print("  L_add_dot = " + L_add_dot ) 

  gsres = True
  gsres@gsMarkerIndex = 1  
  gsres@gsMarkerColor = icla ;; 26 ; 80  ;"red" 
  gsres@gsMarkerThicknessF = 0.008
  gsres@gsMarkerSizeF = 0.03

  do is = 0, ns-1  

     gsres@gsMarkerColor = 1

     if(ieur(is).eq.0) then 
     if(lat(is).le.60.) then 
        gsres@gsMarkerColor = icla 
     else  
        gsres@gsMarkerColor = icle 
     end if  
     end if  

     if(ichi(is).eq.0) then 
        gsres@gsMarkerColor = iclb 
     end if   

     if(iusa(is).eq.0) then 
        print("=== site: "+is+"  "+prtname(is)) 
        gsres@gsMarkerColor = iclc 
     end if   

     if(isam(is).eq.0) then 
        gsres@gsMarkerColor = icld 
     end if 

     if(kbad(is).eq.1) then 
        ;if(kset(is).eq.1) then 
          ;gsres@gsMarkerSizeF = 0.01
          ;print("site name: "+sitename(is)) 
          ;gsn_polymarker(wks,map,lon(is),lat(is),gsres) 
        ;else
          gsres@gsMarkerSizeF = 0.03
          print("site name: "+sitename(is)) 
          gsn_polymarker(wks,map,lon(is),lat(is),gsres) 
        ;end if 
     end if
  end do

  delete(gsres)
   
  end if 

;-------------------------------------------------------------
; Add other things 
;-------------------------------------------------------------

  L_other = False 


  if (L_other) then 

  print("  L_other = " + L_other ) 

     gsres = True
     
     gsres@gsMarkerIndex = NhlNewMarker(wks, "u", 34, 0.0, 0.03, 1.25, 1.0, 0.0) 
     
     gsres@gsMarkerColor = "black" 
     gsres@gsMarkerThicknessF = 0.1
     gsres@gsMarkerSizeF = 0.007
     
     do is = 19, dimsizes(lon) - 1
        gsn_polymarker(wks,map,lon(20:),lat(20:),gsres)
     end do
  
  end if 

  frame(wks)      

  system("ps2pdf "+ploteps+" "+plotpdf) 
  system("mpack -s "+ploteps+" -a "+plotpdf+" k.zhang.iap@gmail.com")  

